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Abstract 

We construct a diffusive matrix model for tfie /?-Wishart (or Laguerre) ensemble for 
general continuous /3 € [0, 2], which preserves invariance under the orthogonal/unitary 
group transformation. Scaling the Dyson index (3 with the largest size M of the data 
matrix as P — 2c /M (with c a fixed positive constant), we obtain a family of spectral 
densities parametrized by c. As c is varied, this density interpolates continuously 
between the Marcenko-Pastur (c — >■ oo limit) and the Gamma law (c — > limit). 
Analyzing the full Stieltjes transform (resolvent) equation, we obtain as a byproduct 
the correction to the Marcenko-Pastur density in the bulk up to order 1/M for all /3 
and up to order for the particular cases /? = 1, 2. 

1 Introduction 

The theory of matrices with random entries, originally devised as a tool to understand 
and predict the spectra of heavy nuclei for which a detailed account of the interactions 
between particles is too complicated, has seen a spectacular resurgence of interest in recent 
years, with a number of unexpected and often surprising applications (see [H El El 15 E] 
for a recent overview). While Wigner and Dyson are usually regarded as the pioneers in 
the field, John Wishart had already introduced random matrices in 1928 in his studies 
of multivariate populations [6]. The Wigner-Dyson (Gaussian) and Wishart ensembles 
(together with a few others) lie at the core of the classical world of invariant matrices, 
characterized by the following main features: 

1. The joint probability distribution (jpd) of matrix entries, collectively denoted by 
P[X], remains unaltered if one performs a similarity transformation X — t- UXU~^, 
with U being orthogonal (real symmetric X), unitary (complex hermitian X) or 
symplectic (quaternion self-dual X) matrix. As a consequence, the eigenvectors of 
such matrices are Haar (uniform) distributed in their respective groups. 



2. The joint distribution of the A'" real eigenvalues P(Ai, . . . , Aiv) can be generically 
written in the Gibbs-Boltzmann form, 

P(Ai, . . . , A^) = ^ exp {-U{Xi, Xn)) (1.1) 
with the Hamiltonian 'H(Ai, . . . , Xn) given by: 

N 

mXi, ...,Xn) = Y^ V{Xi) - /3 ^ In I A, - A^j (1.2) 

i=l j<k 

where V{x) a confining potential that depends on the precise form of the joint distri- 
bution of matrix entries -P[X]. For example, if the entries of X are independent, the 
only allowed potential is quadratic V{x) = /3x^/2, which correspond to the Gaus- 
sian ensembles. If correlations among the entries are allowed, then different potentials 
(all corresponding to rotationally invariant weights) are possible. For example, in the 
Wishart case, V{x) = oo for x < (so that all the eigenvalues are non-negative) and 
V{x) = x/(2(T^) — alogx for x > 0. 

The normalization constant is called the partition function and is simply given 
by the multiple integral 



Zn = J ■■■ JUdX^exp i-niXi, Xn)) ■ (1.3) 



From (1.1), one easily deduces that the system of N eigenvalues of a classically 
invariant ensemble behaves as a thermodynamic system of charged particles arranged 
on the real line, in equilibrium at inverse temperature /3 under competing interactions 



(the confining potential V{x) and the logarithmic all-to-all repulsion term) in (1.2). 
In contrast with the usual canonical ensemble in statistical mechanics, however, the 
so-called Dyson index f3 is quantized and can only assume the values /? = 1,2,4 for 
real symmetric, complex hermitian and quaternion self-dual matrices respectively. 



Lifting the quantization of /3 (Dyson's threefold way) has been a major theoretical 
challenge in view of possible applications e.g. to the quantum Hall effect. Dumitriu and 
Edelman [7j were eventually able to construct ensembles of tridiagonal matrices with 
independent entries whose eigenvalues are distributed as (1.1) with general ^ > 0. Their 
ensemble is however not invariant under similarity transformations, and the eigenvectors 
are not Haar distributed in the appropriate symmetry group. After an earlier attempt in 
the case of 2 x 2 matrices [8] , the explicit construction of an ensemble oi N x N matrices 
displaying at once rotational invariance and a continuous /3 was put forward in [9l [10] 
for the Gaussian ensemble. It was further shown in [9J that only by letting the Dyson 
index /? of that ensemble scale with the matrix size in an appropriate way (namely 
/3 = c/N) one obtains a continuous family of deformed spectral densities parametrized by 
c, interpolating between Wigner's semicircle (typical for /3 ~ 0{1) invariant ensembles) 
and a Gaussian law (properly describing the non- interacting limit /? — )• 0). This result can 
be established in two alternative ways: 
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1. Starting from the stationary joint distribution of eigenvalues (eq. ( 1.1 )), setting /3 = 
c/N ^ and then finding the average density of eigenvalues p{X) = (l/-/V)(X]j 5{\ — \i)). 
In the limit of large A^, this average density can be obtained by a saddle point analysis 



of the partition function Eq. (1.3) in a standard way. Usually, when /? ~ C(l), only 
the energy term ~ 0{N'^) dominates and the entropy term ~ 0{N) is subleading. 
However, when /3 ~ c/iV, both the energy and the entropy terms are of the same order 
(~ 0{N)), which leads to a nontrivial modification of the density. The complicated 
nonlinear singular integro-differential equation for the saddle density, reduces very 
nicely to a Riccati equation for the Stieltjes transform of the density, which can then 
subsequently solved exactly. Finally one obtains the density by taking imaginary 



part of the Stieltjes transform (see section 3.3 for a detailed discussion). 



2. Starting from the dynamical equation of motion of the eigenvalues, one first derives 
directly the equation of motion of the Stieltjes transform of the density via Ito's 
calculus, finds the stationary solution and then obtains the average density by taking 



the imaginary part of the Stieltjes transform (see detailed discussion in section 3.4). 



In section 3.3 and 3.4, in the context of our model, we show that both methods lead 



to the same solution. 

The purpose of this paper is threefold. We first explicitly construct a random matrix 
model l) which is invariant under similarity transformations (and thus has Haar distributed 
eigenvectors), and u) whose jpd of eigenvalues is exactly given by the /3-Wishart ensemble 
of random matrices with a continuous /3 > 0. Then, letting the Dyson index /3 of the 
ensemble scale inversely with the size of the matrix, we analytically derive the density of 
states for this crossover model, written in terms of the Whittaker hypergeometric function 



(see eq. (3.49)), and we show that it continuously interpolates between the Marcenko- 



Pastur law and a certain type of Gamma distribution (see subsection 3.3). Some other 
types of deformations of the Marcenko-Pastur distribution for Wishart-like matrix models 
were reported in the literature (see e.g. |lll I12j). Finally, keeping the Dyson index f3 
unsealed (i.e. it remains of ~ 0{1) for large matrix size A^) but nevertheless continuous 
{(3 > 0), we analyze the full Stieltjes transform equation and we can compute the 1/N 
correction to the Marcenko-Pastur asymptotic density for the /3-Wishart ensemble for all 
values of /3 > 0. Furthermore, using results obtained by Pastur and Lytova in [T^ on the 
noise in the Marcenko-Pastur law, we are also able to derive the correction term in 

the particular cases (3 = 1 and 2. 

The rest of the paper is outlined as follows. In section [2] we introduce the main features 
of the classical Wishart ensemble along with the evolution law for the eigenvalue process. 
In section [sj we construct a 3-parameters matrix model (and the respective evolution 
law for the eigenvalues) that at large times interpolates between the Wishart ensemble 
and so-called CIR processes whose stationary pdf is a certain Gamma distribution. The 
corresponding parametrical density of states is computed exactly in the two ways described 



above (from the saddle point route on the partition function in section 3.3 and from Ito's 



calculus in section 3.4 ) and constitutes a continuous deformation of the Marcenko-Pastur 
distribution (see below). In section [4] we compute the systematic 1/A^ (for all /3 > 0) and 
(for (3 = 1,2) corrections to the Marcenko-Pastur law for the smoothed density in 
the bulk for the (scaled) /3- Wishart ensemble. For the special cases /3 = 1, 2 and 4, the 
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0{1/N) correction term was computed in [T31[TS]. Our result generalizes this to arbitrary 
(3 for the 0{1/N) term and in addition, we obtain the 0{1/N'^) correction for (3 = 1 and 
(5 = 2. We conclude with a summary and discussion in section [5} 



2 Wishart ensembles 

2.1 Real and complex Wishart ensembles 

Let X be a real (respectively complex) Gaussian random matrix of size M x N, i.e. a 
random matrix chosen in the space of M x real (resp. complex) matrices according to 
the law: 

P(X)dX oc exp (^-^Tr(Xtx)^ dX , (2.1) 

where X^ is the Hermitian conjugate of X. In the following, we will denote the real (resp. 
complex) Wishart ensemble by W"^ with /? = 1 in the real case (resp. /3 = 2 in the complex 
case) . 

The real (resp. complex) Wishart Ensemble is the ensemble of (N x N) square matrices 
of the product form W := X^X where X is a real (resp. complex) Gaussian random matrix 
of size N X M. They have appeared in many different applications such as communication 
technology [16], nuclear physics [17], quantum chromodynamics [18], statistical physics of 
directed polymers in random media [T9| and non intersecting Brownian motions \20\ [21} [22] , 
as well as Principal Component Analysis of large datasets |23l |21] • 

The spectral properties of the Wishart matrices have been studied extensively and it 
is known [25] that for M ^ N, all positive eigenvalues of W are distributed via the 
joint probability density function (pdf) 

F,(A..^.,A„)= i.-ii>i:.".^.nA5'"-"«'-'niA.-A/ (2,2) 

i=l i<j 

where Z is a constant normalization factor and where /3 = 1 in the real case (resp. (3 = 2 



in the complex case). Note that the joint distribution defined in (2.2) is in fact well 
behaved for every (3 > 0. 

Another classical result of Random Matrix Theory concerns the asymptotic density of 
states (or spectral measure) for the eigenvalues (Ai, A2, . . . , Aat) of a real Wishart matrix 
W G yV^ in the limit of large matrices, i.e. when N, M ^ 00 with N/M = g G (0; 1] where 
g is a fixed parameter. Let us recall that the density of states of the matrix W is simply 
the probability measure defined as 



1 ^ 



where (3 is introduced for later convenience (/3 = 1 in the present case) and where 
(Ai,...,AAr) are the eigenvalues of W. Setting A = MA, the Marcenko-Pastur Theo- 
rem (see [26]) states that, in the limit N, M ^ 00 with N/M = g G (0; 1], the spectral 
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measure of a Wishart matrix W E converges to a continuous probability density (with 
compact support) given by 



P/3(A) 



1 ^(6-A)(A-a) 



a < A < 6 



27rc72/3g A 
where the edges o, b of the spectrum are given by 



(2.4) 



with again /3 = 1. 

For general /3 > 0, the probability measure is defined again as in (2.3) where this 
time the vector (Ai, . . . , Atv) is distributed according to the law Pg in (2.2). The Marcenko- 
Pastur theorem remains in fact valid for all /3 > in the sense that the probability law 
converges when N, M ^ oo with N/M = g G (0; 1] to the continuous probability density 



Pis in (2.4) for every /3 > 0. 

The probability measure p^ will sometimes be referred to as the spectral density as it 
corresponds to the spectral density of random matrices W G at least when /3 = 1 or 
2. 



2.2 Continuous processes for real and complex Wishart ensembles 

We wish to define here a diffusive matrix process depending on a fictitious time t ^ that 
will converge to the Wishart Ensembles in the limit of large time. The idea is simply to 
set 

Wt := x|Xt (2.5) 

where is a real (resp. complex) random matrix process (of size M x A^) following the 
Ornstein-Uhlenbeck law, 

dXt = —-X.tdt + a d'Bt 

where Bt is a real Brownian (resp. complex) random matrix, i.e. a matrix whose entries 
are given by independent standard Brownian motions. By a standard Brownian motion, 
one means a centered (zero- mean) Gaussian process with covariance function {BtBf) = 
mm{t, t'). 

It is well known that the stationary law of a Ornstein-Uhlenbeck process is the Gaussian 
law and therefore, the real (resp. complex) matrix process Xf converges in law when t — )• oo 
to the law of a Gaussian real (resp. complex) random matrix. Hence, we deduce that the 
real (resp. complex) matrix process defines a diffusive matrix process that converges 
in law to W G with /3 = 1 or /3 = 2. 

It is also easy to check that the positive definite matrix process verifies the following 
stochastic differential equation studied by Bru [27]: 

dWt = -Wtdt + a dBt + a dBl ^/Wt + 1 dt (2.6) 

where is a real (resp. complex) Brownian random matrix and with /3 = 1 in the real 
(resp. /3 = 2 for complex) case. 
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The evolution of the eigenvalue process Xi(t) ^ X2{t) . . . ^ XN{t) is also easy to derive 
8] using perturbation theory to second order 

dXi = -Xidt + 2ay^idbi + (7^/3 ( M + V ^i±^ | dt (2.7) 

V Ai - Afc y 

where 6i's are independent standard Brownian motions and with /3 = 1 in the real (resp. 
/3 = 2 for complex) case. The stationary distribution of the process (Ai, . . . , A7v)(t) is 
necessarily the jpdf Pg(Ai, . . . , A^v) defined in (2.2) (this is true for any /? > and can 
also be recovered using the Fokker-Planck equation for the multivariate diffusion (2.7)). 



3 Crossover between Wishart and CIR processes 

Following [21 [TOj, we aim at defining a diffusive matrix process Wj which converge in 
the limit of large time to a general /3- Wishart matrix, i.e. a matrix whose eigenvalues 



are distributed according to Pjs in (2.2) for general /3 > and with Haar distributed 
eigenvectors. In this paper, we will restrict ourselves to the description of the eigenvalues 
process but the interested reader can find a study of the eigenvectors for a related model 
in ^D|. To simplify notations, we will take in this section a = 1. 



3.1 Preliminary definition: CIR diflFusion process 

We first need to introduce a family of real diffusion processes. Let 5 > be a fixed 
parameter. The CIR process (named after its creators Cox, Ingersoll, and Ross |29j and 
widely used in finance to model short term interest rate) is a diffusion process x{t) defined 
by x{0) := xq > and for t ^ by 

dx(t) = -x{t) dt + 2 y/x(fj dbt + 5 dt . (3.1) 

Using the assumption 5 > 0, it is easy to see that the process x{t) will remain non 
negative for all times t ^ 0. It is also easy to verify that the stationary pdf of the Langevin 



equation (3.1) is the Gamma distribution with shape and scale parameters k = 6/2 and 



9 = 2 defined as ^ 

Ps{x) = — -j;5-ie"t . (3.2) 

2-2T{l) 

In analogy with squared Bessel processes, the parameter 5 will be called the dimension of 
the process x(t). 



3.2 Diff"usive matrix process for general /3-Wishart matrices 

Following [HlllOj, our goal is to construct a diffusive matrix process whose eigenvalues pro- 
cess is asymptotically distributed according to P/j for general f3 E [0, 1]. This construction 
can be extended respectively for general [3 G [0, 2] (resp. /3 G [0, 4]) by using complex (resp. 
symplectic) Brownian motions instead of real Brownian motions in the following. 

We will in fact describe how to handle the value /3 G [0, 1] by using real Brownian 
matrix. This construction can be extended for the values /3 G [0, 2] using complex Brownian 
matrix and also (3 G [0, 4] using symplectic Brownian matrix. 
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The idea is to slice the time interval into small chops of length 1/n and for each interval 
[fc/n; {k + l)/n], to choose independently Bernoulli random variables e'^,k £ N such that 
P[e^ = 1] = p = 1 — lP[e^ = 0]. Then, setting e" = ej^^j, our diffusive matrix process evolves 
as: 

dW]^ = -Wi" dt + dA^ (3.3) 

where the increment matrix c?A" now depends on the value of the additional random 
process e": 

• if e" = 1, then 

dA^ = y^Wf dBt + dBl y^Wf + M Idt. 

where dB^ is an x real[^Brownian increment matrix whose entries have variance 
dt. 

• if ej" = 0, then 

dA]; = dYt + dYj + 5 Idt. 

with 5 > and where dYt is a symmetric matrix that is co-diagonalizable with W" 
(i.e. the two matrix have the same eigenvectors) but with a spectrum given by N 
independent real Brownian increments of variance dt. 

An algorithmic description of how to build (approximatively on a discrete grid) the 
matrix process W" can be found in Appendix [B| 

It is clear that the eigenvalues of the matrix W" will cross at some points but only 
in intervals [k/n; (k + l)/n] for which = (in the other intervals where they follow the 
SDE ( |2.7| ) with parameter /3 = 1 , it is well known that the repulsion is too strong and thus 
collisions are avoided). In this case, the eigenvalues are re- numbered at time t = {k + l)/n 
in increasing order. With this procedure, when ordered A"(t) ^ ... ^ A^(t), we can again 
check as in [9l [10], using perturbation theory, that the eigenvalues will remain always 
non-negative and will verify the Stochastic Differential System (SDS): 

/ \n I \n\ 

dXf = -Xyt + 2^dh + e^M + (1 - e^)6 + E T^T^T^ dt (3.4) 



where the bi are independent standard Brownian motions, which are also independent of 
the process e". 

Note that when e" = 0, the particles A" are evolving as independent CIR processes 
of dimension 5 > as defined in paragraph |3.1[ Therefore, the particles can cross in 
those time intervals, breaking the increasing order so that they will be re-ordered at time 
{\nt\ + l)/n but they will remain non- negative as the dimension 5 is strictly positive. 



Therefore the SDS (3.4) remains well defined at all times t ^ 0. 

One can follow the proof of [10] to prove that the scaling limit (i.e. the limiting process 
when n — )■ oo) of the process (A"(t) ^ . . . ^ satisfies the following SDS 

d\i = -\idt + 2^/\,dhi+ [pM + {l-p)5 + pS^ \'^\'' I dt. (3.5) 



^Here one can use use complex Brownian motions instead to extend the interval of j3 to [0, 2]. 
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One can deduce from the above equation (3.5) the Fokker-Planck equation for the joint 



density P{{Xi},t), for which the stationary jpdf is readily found to be [see the derivation 
in appendix [A] 



P*iXu . . . , A.) = ^e-^ n Af n - (3-6) 



i<j 



The probabihty introduced in (2.2) is recovered here by taking the values p = (3 and 
5 = 0. The corresponding large N, M-limit spectral probability density is therefore given 
by the Marcenko-Pastur law in the case where p = (3 > independent of M. Note that 
with the above normalizations, the spectrum is spread over a region of M+ of width of 
order pM = /3AI. On the other hand, if p = 0, the large A^, M-limit of the spectral density 
is the Gamma distribution with shape and scale parameters k = 5/2 and 9 = 2 (recall 
that it is the stationary pdf of the CIR process of dimension 6): 



1 xl-l. 



A 



Po(A)dA=^;^^A.-^e-^dA. (3.7) 

It is quite natural to ask whether a crossover regime may be found, interpolating between 
the Marcenko-Pastur density {p independent of M) and the Gamma distribution (p = 0). 
A good candidate for triggering such a transition is clearly a parameter p vanishing with 
M as p = 2c/M where c is a positive fixed constant. We discuss this case in the following 
subsection. 

In the next two subsections we compute the crossover density interpolating between the 
Marcenko-Pastur law and the Gamma distribution with shape parameter 6/ 2. This family 
of probability densities is indexed by the three parameters c (such that p = 2c/ M), q = 
N/M and 5^0. More precisely, we compute the limiting density of the probability measure 
Pn = jj YliLi '^Aii when A, M — )• oo with N/M = g G (0; 1] and where (Ai, . . . , Xn) is 



distributed according to the law P* defined in (3.6) with p = /3 = 2c/ M. As mentioned 



in the introduction, the crossover density can be computed via two alternative methods: 



(1) by a saddle point method as shown in subsection 3.3 and (2) by analyzing directly the 



stochastic differential systems introduced in 3.2 above following the analogous route for 



the Gaussian case in Ref. [9]. This is done in section 3.4 We will see that both methods 
yield identical result. 

3.3 Crossover for the spectral density via the saddle point method 



Our starting point in the joint probability law of eigenvalues in Eq. (3.6), where the 
normalization constant (partition function) Z is given by the A^-fold integral 

ndA,e-^"^'>l (3.8) 



'[0,oo]JV ■ ^ 

where the energy function i?[{Ai}] is given by 
E[{^^}] = ^ A, - (|(M - A + 1 - 5) - (1 - 6/2)) ^ In A, - I ^ In I A, - A, I . (3.9) 
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Written in this form, Eq. (3.8) is the Gibbs-Boltzmann canonical weight of a system of 
charged particles on the positive half-line in equilibrium at inverse temperature /? = 1 
under the effect of competing interactions. We wish to compute the average density of 
states {jj: "^^^i 6{X — Aj)) where (O) denote the expectation value of O with respect to the 
probability distribution P* in ( |3.6| ) . There are many ways to compute this average density, 
but the one rather convenient for large is the saddle point route. This was originally 
done by Dyson |30j for the Gaussian random matrices and a physically more transparent 
derivation can be found in Ref. |31j . 

The main idea behind the saddle point method is as follows. In the large N limit, 
the most dominant contribution to the partition function emerges indeed from a set of 
configurations of Aj's that correspond to a particular density p*{X). Naturally then, the 
average computed over the ensemble of Aj's, in this large N limit, will also be given by 
the saddle point density (-^ YliLi ^i'^ ~ -^«)) ~ P*W- It thus suffices to analyze just the 
partition function Z in the large limit and find, in particular, the saddle point density 
p*{X) that maximizes the partition function Z for large N. 

To analyze Z in the large limit, one first defines a 'local' smooth density function 

1 ^ 

'^(^) = ]^E^(^-^*) (3.10) 

i=l 

which is normalized to unity. The main idea then is to split the multiple integration in 



Eq. (3.8) in two parts: First fix the local density p{X) and sum over all microstates (i.e.. 



configurations of Aj's consistent with the local density defined in (3.10)) and then, sum 
(functionally) over all possible local density functions. Roughly speaking, the first step 
corresponds to a partial tracing over microstates by fixing the local density. Notationally, 
on can express this by 

z = [ v[p] [ n^^^^"''"'"^'^['^(^)'{^^}] (3.11) 

where T>[p] denotes a functional integration over the function space and /[p(A),{Ai}] is 
an indicator function that is 1 if the microstate {Aj} is compatible with a given p{X), 



normalized to unity, as defined in (3.10) and otherwise / = 0. The energy function i?[{Aj}] 
associated with a microstate can then be expressed in terms of the local density p{X) using 
the identity J2i fi^i) = ^ I dXf{X)p{X) and one gets ^ 



E[pW] = y / dXXpiX) 



N / dXp{X) In A 



- |iV2 j j dXdX'p{X)p{X') In |A - A'l + |iV j dXp{X) In ^'^^i^j '^^P^^'> ~ ^) 

(3.12) 

where the last term includes a Lagrange multiplier Ci that enforces the normalization of 
the local density to 1. The next-to-last term accounts for the self-energy term (A — t- A') 
that needs to be subtracted. Note indeed that in the original discrete sum Yli^j 
the eigenvalues do not coincide. This means that the integral over A and A' should exclude 
the region where |A — A'| is less than the typical spacing between eigenvalues, which is 
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proportional to \/N p[X) |30j . The contribution of this thin sHver is the next-to-last term, 
up to an additional contribution that can be absorbed into Ci. 



Once this is done, equation (3.11) simplifies further and one gets 



j P[p]e-^W^)1 J[p{\)] 



(3.13) 



where J[p{X)\ = oo]^ Ilj '^'^i-^ {-^«}] an entropic factor that just counts how 
many microstates are compatible with a given local density function p(A). This can be 
estimated very simply by the following combinatorial argument involving arrangement of 
N particles in K boxes. Let us first divide our one dimensional line into K small boxes 
of equal width. We have N particles that need to be distributed into the K boxes with 
occupation numbers {^1,71,2, . . . , nx}- The number of ways this can be done is simply 



ni\n2 \ ■ ■ ■ riK^. 



(3.14) 



Setting Pi = rii/N (the local density in box i) and using Stirling's approximation A^! 
j\^N+i/2 ^-N ^^g^j-^g ^]^g fact that ^^rii = N), we have: 



Nl 



ni!n2! • • • n^^! 



~ e 



(3.15) 



which in the continuum limit, becomes ~ f '^^pW^^pW ^ Thus, the entropic factor can 
also be expressed as a simple functional of the local density p{X). Inserting this expression 
in the functional integral over the density, yields: 



(3.16) 



where the free energy F[p(A)] is given by: 



P-N 
2 



f((l-l)Ar + l-5 



dXdX'p{X)p{X') In |A - A'l + ( 1 - I 



-0 ldXp{X)lnX 
dXp{X)lnp{X) + Ci 



dXp{X) - 1 
(3.17) 



Note that for p ~ 0{1/N) the entropy term becomes of the same order of the energy 
term, while in the usual case p ~ 0{1) the entropy contribution is subdominant in the 
large N limit and is therefore disregarded. 

Setting now p = 2c/M = 2cq/N, we get: 



F[p{X)] 



1 



dXXp{X) 



cq 



1 



cq 



dXdX'p{X)p{X')ln\X- X'\ + (l 



cq 

N 



dXp{X) In X 
dXp{X)lnp{X) + Ci 



dXp{X) - 1 
{3.1. 
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We set a = cq{l/q — l) — (1 — 5/2) and take A'' — )• oo (so that the term cq/N drops out). 
Now, the saddle point density p*(A) is obtained by minimizing the free energy F[p[X), i.e., 
by taking the functional derivative ^ = and p{\) == /o*(A). This gives the saddle point 
equation 

^-a\u\-2cq j dA'p* (A') In |A - A'| + Inp* + = (3.19) 

where C2 = Ci + 1 is just a constant. For notational simplicity, in the rest of the subsection 
we will denote the saddle point density p*{\) simply by p{X). 



Taking one more derivative of (3.19), we get 

P(A') 



1 a 

- - - - 2cg Pr 

2 A ^ 



A- A' 



(iA' + 



Next, we define the Stieltjes transform: 

H{z) = 



dX 



for z complex and outside the support of p. By definition, for large \z\, H[z) 

integrating over 

p{x)dx f p{x': 



Multiplying eq. (3.20) by p{X)/{X — z) and integrating over A, we have: 

p{X)dX 



P(A) 
A- 



dX — a 



X{X-z) 



2cq Pr 



X-z 



A- A' 



■dX' + 



p'{X)dX 
X-z 



(3.20) 



(3.21) 
-l/z. 



(3.22) 



and we analyze each of the four contributions separately. 



1. Ti 

2. T2 



1 J 



X-z 



dX 



\H{z) 



-a J x{\^-z) ■ rewrite this as: 



a / p{X)dX 



X-z 



implying: 



where bi = a J dA^^ 



T2 = — H{z) H 



(3.23) 
(3.24) 



3. T3 = —2cq Pr J ^'^^^^ f dX' which we rewrite as: 



T3 = 2cq 



dXdX'p{X)p{X') 



1 



1 



= 2cqH'^{z) - 2cq Pr j 
By renaming A — t- A' and A' — t- A, we get: 

T3 = 2cqH^{z) + 2cq¥T j 
= 2cqH\z)-T3 

Solving for T3 we get: 



X-z A - A' 
dXdX'p{X)p{X') 
(A-A')(A'-z) 



A' 



dXdX'p{X)p{X') 
{X-X'){X-z) 



n = cqH'{z) 



(3.25) 
(3.26) 

(3.27) 
(3.28) 

(3.29) 
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p'(X)dX 



which we integrate by parts, obtaining: 



1 



X-z 

Cl 



+ 



P(A) 

{x-zy 



:dX 



z 



+ H'{z) 



(3.30) 



In the derivation above, we assumed 61 and ci to be finite. This is not completely 
obvious, because p{X) at an edge point may diverge. However, by imposing that for large 
z, H{z) — )• —1/z, it is immediate to derive that bi + a = 1/2. Thus, one may regularize 
the density near the edge points so that 61 and ci exist individually, but eventually their 
sum is universally 1/2 and hence is independent of the specific regularization near the 
edge. 

Adding up the four contributions, we get the equation: 

dz z z 2 

Thus we find the following differential equation for the Stieltjes transform H: 



dH ^ o 1 / a\ 1 



(3.31) 



where we have set 



a = (2 - (5) - 2c(l - g), 7 = eg. 



In the next subsection, we will derive the same equation via Ito's stochastic calculus route. 
The density p{X) (normalized to unity) can then be read off from 



p{X) = -\m[H{z A)] 
vr 

where z — t- A occurs inside the cut on the real axis. 



(3.32) 



To solve the Riccati equation (3.31), we make a standard substitution 

H[z) = -'^ = -dAnu{z). 
7 u{z) 7 

This gives a second order differential equation for u{z) 



u(z) + 



1 



1 + 



a 



u'{z) + ^u{z)=Q. 



It follows from Eq. (3.33) and the asymptotic behavior of H{z) that 



u{z) 



|2|-s>oo 



(3.33) 



(3.34) 



(3.35) 



where is a constant. 
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tion 



To reduce Eq. (3.34) to a Schrodinger like differential equation, we make the substitu- 

(3.36) 



u[z) = e 

and we find the following equation for ip 



1 1 — a a / a 
16 z 8 4 V ~ 4 



1 



0. 



Making further a rescaling 'ip{z) = y{z/2), it reduces to the standard form of the Whittaker 
differential equation ||32j 



y"{z) + 



1 A 

■7 + - + 
4 z 



z^ 



(3.37) 



(3.38) 



where 

a l a 

Note that the solution of this differential equation does not depend on the sign of ^, hence 
we take the absolute value. The differential equation (3.37) has two linearly indepen- 
dent solutions W^^^(z) and VF_^^^(— z). The Whittaker function W(^^^{z) has the following 
asymptotic behavior [32] 



Thus the general solution of u{z), using Eq. (3.36), reads 

l{z) = e-^l^ z^l^ [Ci W^c,m(^/2) + C2 W.^^^{-zl2)\ 



UK 



(3.39) 



(3.40) 



where C\ and C2 are arbitrary constants. Using the asymptotic behavior in Eq. (3.39) it 
is easy to check that only the second solution has the right asymptotic behavior in Eq. 
(3.35). Thus, finally, we have our solution 

u{z) = C2 e-'/' z^l^ W.^,^{-z/2) (3.41) 



where C, and ^ are given in Eq. (3.38). 

By plugging this solution (3.41) into Eq. (3.33) and using Eq. (3.32), we find the 
following expression 



P(A) 



C2 (Im(Ty:^^^)Re(Ty_c,^) -Im(W^_c,^)Re(VF:^_^))(-A/2) 



(3.42) 



2^7 |T^_c,^(-A/2)|2 

where Re and Im denote respectively the real and imaginary parts. Using the linear dif- 
ferential equation verified by the Whittaker functions (3.37), it is easy to see that the 
derivative with respect to A of the Wronskian type function {lm{W'_^ ^)Re(Ty_(^^^) — 
Im(W_^^^)Re(PFi^ ^)) is equal to 0. Hence the Wronskian appearing in the numerator 
in (3.42) is simply a constant. 

Collecting all the constants together, we get: 

A 



P(A) 



|Ty_C,;.(-A/2)| 



2 ■ 



(3.43) 
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The overall normalization constant A has to be fixed from p(A) dX = 1. Thus we get, 
after rescaling A/2 — )• A, 

(3.44) 



1 
A 



lo |W^-0«(-A)p- 

This integral in Eq. (3.44) can be done in closed form. First, we first use the well known 
identity [33] 

W^,^{z) = u{ij-C + 1/2, 1 + 2/i; z) (3.45) 

where U{a, b; z) is the Tricomi hypergeometric function (or Kummer function) that behaves 
for large z as U{z) ~ z~°'. Using this in Eq. (3.44) gives 



1 f'°° 

- = 2 / dXX-^^'-^e-^\U{fi + C + 1/2,1 + 2^i;-X)l 
A Jo 



(3.46) 



It turns out that there exists an interesting integral representation in a paper by Ismail 
and Kelker 



1;^ 



°°dte-H-\^^^ , ,,2 , ,lU(a,b 

\U{a,b;-t)\-^ = T{a)T{a-b+2^ ^' 



z + 1 



z U{a,b;z) 



for a > 0, 1 < 6 < a + 1 
(3.47) 

Note that in Ref. [M] they use the notation ijj(a,b,z) instead of U{a,b;z), but it is the 
same function. Our ^ and C satisfy the condition of validity of this identity: a > and 
1 < 6 < a + 1. Taking z — )■ oo limit on both sides and using U{z) ~ we arrive at the 
following exact expression of the normalization constant 



1 

A 



2r(/i + c + i/2)r(c-// + 3/2) 



(3.48) 



This leads to the following final expression for the spectral density, which is the central 
result of our worlo 



Pc(A) 



1 



2r(/i + c + ^)r(c - Ai + 



\w. 



1 



(3.49) 



with the following values for the parameters 

{2 - 6) - 2c{l - q); ( = cq 



a 



a 



— ; and n = -\a 
4 4 



The above expression is the analogue, in the present context, of the Askey-Wimp-Kerov 
one-parameter family of models found in [9] , that smoothly interpolates between the Gaus- 
sian distribution and Wigner's semi-circle. 

Let us now consider the limiting case c — )• first. In this case, we have a = (2 — 5), 
( = 5/4- 1/2 and fi = 5/4. Thus, VF„^,^(-A/2) = W^i/2-5/4,5/4(-A/2). It turns out 
that for these special values of the indices, the Whittaker function simply reduces to 
^1/2-^/ 4,5/4 ( ~ A/2) oc A"''/^"'"^/^ e"^/^ up to a proportionality constant t32j. Substituting 
this in (3.49), we then recover the CIR density in Eq. (3.2). The limit c — )• oo is more 



■^We add the subscript c in the notation of the density p to recall the dependence in c; the density p 
depends also on the two parameters 5 and q but wc omit to subscript those. 
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tricky as one needs to rescale A — ?• Ac and take the large c limit carefully. This can be 
done and one recovers the Marcenko-Pastur law. This can be rewritten (without rescaling 
A) as 

1 \/(7+ - A)(A-7_) 

Pc(A) -c^.^ — 1{7_<A<7+} > (3-50) 

where 7-1- = 2c (1 it y^)^. 

Using standard results on Whittaker functions (see e.g. [32l[33]), it is easy to compute 
the asymptotic behavior of /9c(A) for A — )• and A — )• +00. Up to multiplicative constants, 
we have 

Pc(A) ~A^+oo e~2 , 

and 

We plotted in Fig. [ijthe density pc for c = 0, 1, 2, 3, 4, 5, 10 and q = 1/2,5 = 1, showing 
the progressive deformation of the Gamma distribution with shape parameter 5/2 towards 



the Marcenko-Pastur distribution (3.50) with parameter q = 1/2. The critical value of c 
at which the divergence at A — )• 0+ changes to convergence is c* = (2 — 5)/ {2(1 — q)) = 1. 
As expected, in Fig. [T| the curve with second highest value at the origin corresponds to 
c = 1 and converges when A — ?■ 0+ to 1/2. The curve with highest value at the origin is 
the Gamma distribution with shape parameter 5/2 = 1/2 and diverges at 0+ to +00. The 
other curves corresponding to c = 2, 3, 4, 5, 10 converge to when A — )• 0+. 



We have also verified our analytical result for the crossover density in (3.49 ) numerically 
for the sample value c = 1 and found very good agreement (see Fig. [2]). 

3.4 Crossover for the spectral density via Ito's stochastic calculus 



In this subsection, we want to re-obtain the result Eq. (3.31) of the previous subsection 



via Ito's calculus. We therefore consider the process (Ai(i), • • • , AAr(t)) which verifies the 



stochastic differential system (3.5 ) with the scaling relation p = 2c/M. The idea is to work 



out the evolution equation of the probability measure 

1 ^ 

i=l 

in the large N limit. We expect the equilibrium of this evolution equation to be the solution 



of (3.31) 



In the following, / is a smooth function. Using Ito's formula for J f(x)p*j^{dx), Eq. 
(3.5) and the scaling relation p = 2c/M, we obtain (see [35] for similar calculations) 

d I f{x)p%idx) = j {-x + 2c+{l- ^) 5^ f'{x)p%{dx)dt (3.52) 



c 



+ 2(1--) / xnx)pUdx)dt 



+ cqff [ /'(^) f'(y^ + y)p%{dx)p%{dy)dt + dM^ , 
Jo J J x-y 
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Figure 1: Density pd^) for c = 0,1,2,3,4,5,10 of Eq. (3.49) showing the progressive 
deformation of the Gamma distribution (3.2) with parameter 5 = 1 towards the Marcenko- 
Pastur distribution with parameter q = 1/2. The value /Jc(0) at the origin decreases when 
c increases. 
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where dM, 



N 



N 



Yli=i f i''^i)dbi is a noise term of variance / \f' {X)^ p^j^{d\)dt. 



When N, M ^ oo with N/M = q, this noise term is of order 1/y N . 

In the large N,M hmit, the stationary probabihty measure p solution of Eq. (3.52) 
therefore satisfies to leading order (keeping only the terms of order 1) 



(3.53) 



{-X + 2C + S) f'{x)p{dx) + 2 j xf"{x)p{dx) 
f f ^'^''l~^'}^\ x + y)p{dx)p{dy) = 



x-y 



Applying Eq. (3.53) to the particular function f{x) = for z G C\M and denoting H{z) 
the Stieltjes transform of the probability measure p, we obtain the following differential 
equation for H 



[H{z) + zH'{z)] - (2c + 6) H'{z) + 2 [2H' {z) + zH" {z)] 
+2cq [H'{z) + H^{z) + 2zH{z)H'{z)] = . 



(3.54) 



Eq. (3.54) can be rearranged as 

2cqH{z)[H{z) + 2zH'{z)] + \[H{z) + 2zH'{z)] + ]^H{z) + [3H'{z) + 2zH"{z)] 

+ [{l-5)-2c{l-q)]H'{z)=Q. (3.55) 

Eq. (3.55 ) can be integrated easily by doing the change of function G{z) = zH{z'^). Indeed, 
we just need to write (3.55) for z^ instead of z and then multiply the corresponding 
equation by z to obtain the following equation 



2cqG{z)G'{z) + ]^[zG'{z) + G{z)] + \g"{z) + 



1 



c{l-q) 



2zH'{z'^) =0, 



which can be integrated as 



2cqG'^{z) + 



z + 



{l-6)-2c{l 



G{z) + G'{z) = -1 



(3.56) 



where the integration constant is chosen so that zG{z) ~ —1 when \z\ — )• oo. Note that 
the asymptotic behavior for H is therefore also zH ~ — 1 when \z\ — )• oo as is expected for 
the Stieltjes transform of a probability measure. Rewriting now Equation (3.56) in term 
of the function H, we obtain exactly Eq. (3.31). 



4 Correction to the Marcenko-Pastur law for large but finite 
dimension 

In this section, we come back to the case of generalized Wishart matrices for which particles 
are distributed according to the jpdf Pj3 with general parameter /? > (not scaling with 
M). We want to compute the first correction terms to the Marcenko-Pastur density p^ 



(defined in Eq. (2.3)) for large but finite N,M with N/M = q e {0; 1]. 
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We are therefore interested in computing functionals of the form (J f{X)p^{d\)) where 
/ is a test function and where (Ai, A2, • • • , Aat) is distributed according to the jpdf 
defined in Eq. (2.2) for (3 > 0. The idea is to use the stochastic process Xi{t) ^ ... ^ XN{t) 
following the SDE (2.7) that converges in law when t — t- 00 to the vector (Ai ^ . . . ^ Aat) 
distributed according to P^. 

Note that in order to have a well behaved limiting spectral distribution with edges that 
do not depend on M or on /3, we will choose in this section a = l/\/Mf3 (or alternatively 
rescale all eigenvalues by l/M/3). 

Using again Ito's formula for f /(A)/9^(dA) (where is still defined by Eq. ( 3.51| )) 
and Eq. (2.7), we obtain 

d [ fix)pUdx) = [i-x + 1) f'ix)p%idx)dt + -^(2-f3) [ xf"{x)p%{dx)dt (4.1) 



+ 



Af/3 

nx)-ny) ^^ ^ y)pt (^dx)p%{dy)dt + dMi" 

x-y 



where dM^ 



N 



7fJ2iLi'^\/ f'i^i)dbi isanoise term of variance ^^^^^ / Xf'{X)'^p%{dX)dt. 



Let us introduce the Stieltjes transform Ht{z) of the probability measure p'j^ defined 



as 



p%idx) 



X — z 



We now apply (4.1) to the particular function /(x) = l/{x—z) and we take the expectation 



with respect to the bf, Eq. (4.1) rewrites as 
d{Ht) 



at 



m + z 



OjHt. 
dz 



+ 



d{Ht 
dz 



+ {Hty + 2z{H-, 



dz +M/3^^ ^> 
d{Ht 



dz 



dz"^ 



dz 



d 



+ q[{H^)-{Htf]+zq— [{H, 



(4.2) 

{Htf] 



where (•) denotes the expectation with respect to the Brownian motions 6j. The two last 
terms come from the replacement of {H^) by (Ht)'^ in the third term of the right hand 
side of (|4.2|). By setting 



Ftiz) = mzf) - {Htiz))\ 



Eq. (4.2) can be rewritten as 

d{Ht) 



dt 



{Ht)+z- 



dz 



dz 



M/3 



+ 



dz 



^ + {Htf + 2z{Ht 



dz 



+ q 



Ft + z 



dz 
dFt 



dz"^ 



(4.3) 



(4.4) 



dz 



To simplify notations, we will now omit the (•) and write H instead of {H). The stationary 

d'^H 



solution of (4.4) writes simply as 

dH' 



H + z 
+q 



dz 
dH 



dH I , 

+ —(2-/3) 



dz 



Ml3 
dH 



^dH 

dz dz'^ 



(4.5) 



^ +H' + 2zH- 
az dz 



+ 



F + z 



dF 

dz 



0. 
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which can be rewritten as 



+ 



1 



2M/3 



(2-/3) + (-l + 



+ 



F + z 



dF_ 

dz 



0. 



Eq. (4.6) can be integrated easily by doing the change of function G{z) = zH{z ). 
Indeed, we just need to write (4.6) for instead of z and then multiply the corresponding 

(4.7) 



equation by z to obtain the following equation 



qG 



dC 



1 

dz^2 
1 



dC ^ 
Z— + G 
dz 



+ 



1 



+ 



1 



2M/3 



4M/3 
2-/3) + (-l + 



(2-/3) 



d^G 
dz^ 



dH , n. 



+ qz 



0. 



Equation (4.7) can be straightforwardly integrated with respect to z as 

G{z 



qG^ + zG + 



1 . n.dG 



2M/3 



(2-/3) + /3(-l + g) 



-1 (4.8) 



where the integration constant is chosen so that zG{z) ~ — 1 when \z\ — )• c«. Note that 
the asymptotic behavior for H is therefore also zH ~ —1 when \z\ — ?• cxd as is expected for 
the Stieltjes transform of a probability measure. Rewriting now Equation (4.8) in term of 
the function H, we obtain 



qH^ + H 



1 + 



1 



1 + 



1 2-/3 
M f3 



I 2-13 dH I , , , , 



Now, using the result about the noise in the Marcenko-Pastur law obtained in [13j by 
Lytova and Pastur, we know that in the limit of large M, we have for /3 = 1 or 2, 



F(z) 



(4.10) 



1 1 1 

Wq^2p^ 



7+ 



7+ 



dXdfx 



iq-{X-il+q))if,-il + q)) 



(A - z)2(^ - z)2 _ (A - (1 + q)yy/4q - (/x - (1 + q)^ 



with 7± = (1 ± V^f- 

The idea to obtain the correction to the Marcenko-Pastur law is to use perturbation 
theory in Eq. (4.9). More precisely, we want to compute explicitly the coefficients po,Pi 
and p2 such that the eigenvalue density of a /3-Wishart matrix writes under the form, in 
the hmit of large A^, M with N/M = q, 



p{X) = poiX) + ^Pi(A) + ^P2(A) + 



(4.11) 



Note that this asymptotic expansion (4.11) is obtained by perturbation theory and 
therefore is valid only for the values of A such that the correction terms pi{X)/M and 
P2{X)/M^ are negligible compared to the leading term po{X) in the limit of large M, i.e. 
for the values of A such that po{X) 7^ 0. The expansion (4.11) is not valid outside the 
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Marcenko-Pastur sea, i.e., it breaks down near the edges (see below). In addition, here we 
are talking about smoothed density, hence it contains no oscillatory term in the finite N 
bulk corrections [14j . 

To this purpose, we first write H{z) under the form 



H{z) = Ho{z) + ^Fi(z) + ^H2{z) + 



(4.12) 



and we plug Eq. (4.12) into Eq. (4.9). By solving the equation to leading order, we find 
the following expression for Hq{z) 

1 -{z + q-l) + v/(^-7-)(^-7+) 



^o(^) 



2q z 
1 -{z + q-l) + ^{z-{l + q)f-Aq 
2q z 



(4.13) 



We deduce from this the famous Marcenko-Pastur result: the eigenvalue density converges 
in the limit of large N,M with N/M = q to the Marcenko-Pastur density as expected 
given by 



Po(A) 



1 V(A-7-)(7+-A) 



2TTq A 

For all value of /3 > 0, we can now compute the 1/M correction to the Marcenko-Pastur 
density by plugging Eq. (4.12) into Eq. (4.9) and solve to order 1/M. This gives the 
following expression for Hi{z) 



Hi{z) 



l\ 1 



1 



+ 



1 



,z-7+ z-7_y ^(z-7+)(z-7_) 
and the corresponding 1/A^ correction to the density is then given (for all /3) by 



Pi(A) 



2) q 



1 



(<5(A-7+)+5(A-7-)) 



1 



vr 



V(A-7-)(7+-A) 



(4.14) 



Comparing Hi{z)/M with Hq{z), we see that the correction term ceases to be negligible 
when |A — 7-1-1 ~ M^^/^, as expected: this is indeed the standard edge scaling that defines 
the Tracy- Widom region. Note that this 0{1/N) correction term was derived earlier [ |14j 
and references therein] for /3 = 1, 2 and 4, but our result is valid for general /3. 

For the particular value /3 = 1 or 2, we can use the result of Lytova and Pastur stated 
above in Eq. (4.10) to compute the 1/M^ correction with the same method by solving the 
equation until order 1/M^. 

Let us first compute an explicit expression for F{z) from the integral representation 
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in Eq. M.IO) 



1 



4g 



2/37r2 



7+ ]_ 



7- (A-^)V4'7-(A-(l + g))2^ 

2- 

+ o{l) 



A - (1 + 



1 

2/37r2 



(A-z)2 ^4g_(A-(i + g))2 
4g 



2z - 7- - 7+ 

{z - 7_)(z - 7+) y/(z-7_)(2;-7+) 



1 (2z - 7- - 7+) 
(z - 7_)(z - 7+) 2 (z - 7_)(2; - 7+) 



2 \ 2 



+ 0(1) 



1 



1 



2/3 

1 1 

(^-7-) (7+-^) 



4g- 



(2z - 7- - 7+)2 

(z-7_)2(2;-7+)2 

+ o(l). 



^ _ 1 (2z - 7- - 7+ 

2(^;-7-)(^-7+) 



2 \ 2' 



+ 



Then we can turn to compute H2{z) and deduce from this computation the expression 
for p2{X) 



P2W = -^{l-l) ^ 



1 



2 /3y <z^(A-7-)(7+-A) 



1/1 1 

+ 



2 VA-7- A -74 



+ 



2 V(A-7-)2 (A -7+) 



+ 



2g/3 (A-7_)3/2(7+ - A)3/2 



Again, the comparison of this correction term with the dominant term indicates that 
our perturbation expansion breaks down when |A — 7-1-I ~ M~2/3^ 



5 Conclusions 

In summary, we proposed a random matrix model (invariant under similarity transfor- 
mations) whose joint density of eigenvalues is given by the classical /3-Wishart ensemble 
where the quantization of the Dyson index /3 is lifted. The procedure is constructive and 



is described in section 3.2, The resulting ensemble is by construction invariant under simi- 
larity transformations with Haar distributed eigenvectors The diffusive evolution equation 
for the eigenvalues involves the Dyson index of the ensemble as a free parameter. Letting 
it scale with the size M of the matrix, the spectral density of the ensemble becomes a 
one-parameter continuous family interpolating between the familiar Marcenko-Pastur dis- 
tribution and a certain type of Gamma distribution. On the other hand, keeping the Dyson 
index unsealed but not quantized, we showed that a careful analysis of the full Stieltjes 
transform equation lead naturally to 1/N and corrections (and possibly systemati- 

cally to any order) to the average spectral density (Marcenko-Pastur) for all /3 7^ 2. This 
then extends the previous work |14^ I15j on the 0{1/N) correction term for /? = 1, 2 and 
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4. To order 0{1/N'^), our result is valid for /3 = 1 and /3 = 2. It would be interesting to 
see (or conjecture) if this formula to 0(1/ N'^) term is valid for general /3. 

In this work, we have computed the crossover density as a function of the interpolating 
parameter c. In the limit c — >• oo, it reduces to the standard Marcenko-Pastur density, 
whereas the opposite limit c — )• corresponds to the Gamma laws associated with the CIR 
process. It would be interesting to extend our analysis to the distribution of the largest 
eigenvalue. As in the case of bulk density, we would expect a c-dependent distribution 
for the largest eigenvalue, properly centered and scaled, interpolating between the Tracy- 
Widom distribution (c — )• oo limit) and Gumbel distribution (as c — )• 0). 



A Derivation of (3.6) 



The Fokker-Planck equation for the transition probability density P{\i, ■ ■ ■ , Aat; t) of the 
process (Ai(t), • • • , AAr(t)) which satisfies the stochastic differential system (|3.5|) reads 



dP 
~dt 



N 

E 

i=l 



_d_ 



P -Xi+pM + {1-5)p + pY^ 



At + A;, 
A,; — Xi 



& 



i=l 



The stationary solution is the solution which does not depend on time t, satisfying 



N 

E 

i=l 



_d_ 
dXi 



P -Xi + pM + {1 - 5)p + pY^ 



kj^i 



Xi + Xk 
Xi — Xk 



i=l * 



It is easy to check using elementary algebra that the jpdf P* defined in (3.6) verifies Eq. 



(A.l) as in fact we can verify that for all i, 



P* -Xi+pM + {l-5)p + pY, 



k=/=i 



Xj + Afc 
Xi — Xk 



(A.3) 



B Algorithmic description of how to build the process W 
in practice 



Let us describe shortly an algorithmic description of how to build the process W" in prac- 
tice, on a discrete grid. First note that this algorithmic description needs a discrete grid 
and that it does not reproduce exactly the process WJ^ but only a discretized approxima- 
tion of it. Choose a large value of n and an initial symmetric matrix Wq. The construction 
is iterative. Suppose that the process is constructed until time k/n and let us explain how 
to compute the matrix W"^_|_^.jy^ at the next discrete time of the grid, {k + l)/n. 



1. Step 1. We first need to compute the matrix yW^^^. It suffices to compute the 
orthogonal matrix O'^^^ such that 



^k/n — ^k/n^k/n*-'k/n 
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where 5]^^^ is the diagonal matrix composed of the eigenvalues of W^^^ (in increasing 
order) . The eigenvalues of the matrix W^^^ should be non negative as the eigenvalues 
process of W" are almost surely non negative at all time t. However, due to the 
discretization scheme necessary for algorithmic procedure, the non-negativity can 



fail. To avoid this problem, we define y^fc/n ^ 



where y/S^y^ is the diagonal matrix composed of the square roots of the absolute 
values of the eigenvalues of W^^^ (again in increasing order). 

2. Step 2. We sample the Bernoulli random variable with P[e^ = 1] = p = 1 — P[e^ = 
0]. 

3. Step 3. It depends on the value of e^: 

• if = 1, we sample a N x N matrix G„ filled with independent Gaussian vari- 
ables with mean and variance 1/n and then we compute the matrix W"^^^^y^ 
by the formula 

Wr.+i)/„ = (l - ^) W^/n + ^/W^ G„ + Gt + Im I . 

• if = 0, we sample N independent Gaussian variables {zi, ■ ■ ■ , zn) with mean 
and variance 1/n. We then compute the matrix Y„, which is co diagonalizable 
with the matrix W^^^, defined as the product 

Yn := 0^/„Diag (zi, Z2, . . . , zjv) O^/J . (b.2) 
Finally we obtain the matrix W"^^^^^ by 
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